Deriving Image Data from Data Received from an Ultrasound Probe

ABSTRACT

Apparatus for producing image data representing a specimen being imaged using ultrasound includes a processor and memory. The processor is configured to receive input data representing the output from an ultrasound probe, process the input data to produce output image data, and output the output image data. The processor is configured to carry out the step of producing output image data by obtaining a first matrix that is dependent upon attributes of the probe and on the way in which the input data was produced, initialising estimates for the image data, the echogenicity of the image data, and the variance of the noise of the image data, and iteratively performing the following steps until a predetermined level of convergence is reached: (a) updating the estimate of the image data using: the first matrix, the estimate of the echogenicity, the estimate of the variance of the noise, and the input data, (b) updating the estimate of the echogenicity using the estimate of the image data, and (c) updating the estimate of the variance of the noise using the estimate of the image data. Step (c) is performed using an update rule that assumes the noise to be varying across the image data but locally invariant for any small enough region of the image data. The processor then uses the converged estimate of the image data to produce the output image data.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to an apparatus and method for producingimage data representing a specimen being imaged using ultrasound.

2. Description of the Related Art

Ultrasound machines produce an image of a specimen, usually an areabelow a patient's skin, by using a probe including one or more acoustictransducers to send pulses of sound into the specimen. The pulses arereflected back from materials within the specimen and the receivedsignals are used to generate an image that can be displayed to anoperator and/or stored.

However, the pulses generated by the transducers have a focal depthoutside which the images become increasingly blurred. Thus, it isdifficult to accurately image a specimen that has a large depth.

BRIEF SUMMARY OF THE INVENTION

According to an aspect of the present invention, there is providedapparatus for producing image data representing a specimen being imagedusing ultrasound as set out in claim 1.

According to another aspect of the present invention, there is provideda method of producing image data representing a specimen being imagedusing ultrasound as set out in claim 11.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows an ultrasound machine suitable for use with the invention;

FIG. 2 shows components of the ultrasound machine shown in FIG. 1;

FIG. 3 illustrates a probe shown in FIG. 1;

FIG. 4 illustrates an array of transducers shown in FIG. 3;

FIG. 5 illustrates a beam formed by a plurality of transducers shown inFIG. 4;

FIG. 6 illustrates how the plurality of transducers shown in FIG. 5create the beam;

FIG. 7 details steps carried out by the ultrasound machine shown in FIG.1;

FIG. 8 is a diagram of the process of producing image data by theultrasound machine shown in FIG. 1;

FIG. 9 details steps carried out during FIG. 7 to produce image data;

FIG. 10 illustrates the space being imaged;

FIG. 11 is a diagram of a pulse being reflected from a scatterer;

FIG. 12 is a diagram of a plurality of pulses being reflected from ascatterer;

FIG. 13 shows equations used during the steps shown in FIG. 9 to obtainthe image data;

FIG. 14 shows equations used during the steps shown in FIG. 9 to obtaina first estimate of the image data;

FIG. 15 illustrates a transformation into a different domain tocalculate a matrix shown in FIG. 14;

FIG. 16 is a diagram of the algorithm used to obtain the image data;

FIG. 17 details steps carried out during FIG. 9 to carry out thealgorithm shown in FIG. 16;

FIG. 18 shows equations used during the steps shown in FIG. 17 to updatean estimate of the image data;

FIG. 19 details steps carried out during FIG. 17 to update an estimateof the image data;

FIG. 20 shows equations used during the steps shown in FIG. 17 to updatethe estimate of the covariance of the noise of the image data; and

FIG. 21 illustrates an alternative ultrasound machine that may be usedwith the invention.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS FIG. 1

FIG. 1 illustrates a typical ultrasound machine. This includes aprocessing system 101 that receives input from a plurality of probes102, 103, 104 and 105. A probe is held against a specimen to be imaged.Usually, this is the skin of a patient and the area immediately beneaththe skin is imaged, but ultrasounds may also be taken of organs thathave been removed from the body or of other materials.

A pulse is generated by the processing system 101 and sent to the probein use, for example probe 102. This pulse is sent through the specimenand is reflected back. Data representing the reflections is returned toprocessing system 101 and displayed on display 106. An interface 107 isprovided to allow an operative to control the machine.

Many alternative forms of ultrasound machine are available. Inparticular, portable machines that resemble laptops are often used, butthese generally have less processing power and less memory, which canlimit the quality of the displayed images.

FIG. 2

A block diagram of the processing system 101 shown in FIG. 1 is laid outin FIG. 2. A control unit 201 controls the overall process. It instructspulse generator 202 to generate the appropriate pulse which is sent toamplifier 203 before being forwarded to a switch 204. The pulse is thensent to a probe, such as probe 102. When data is received back from theprobe, this is sent via switch 204 to a time gain control unit 205.Analogue-to-digital conversion of the data is performed in unit 206 anda processing unit 207 converts the data into image data that is sent todisplay unit 106. Control unit 201 is then informed to start the processof generating a pulse again.

The processing system shown in FIG. 2 differs from traditionalultrasound processing systems, which generally perform envelopedetection before the analogue-to-digital conversion, in that there areno non-linear steps before output is displayed. This reduces the amountof memory required.

Other processing systems that include a unit or units capable ofcarrying out the steps of the invention may also be used.

FIG. 3

Probe 102 is represented in FIG. 3. It is connected to processing system101 by a cable 301 down which a control wire 302 is passed. This controlwire sends signals to and receives signals from an array 303 ofpiezoelectric transducers.

A coating 304 forms the outer lower surface of probe 102. This canassist with impedance matching. The probe is placed against a patient'sskin 305, usually with a layer of water-based gel 306 between thecoating 304 and the skin 305.

When the probe receives signals from processing system 101, it producessound waves using transducer array 303 as instructed, and these arepassed into the patient's body as shown by arrow 307. Echoes 308 arereceived back to the transducers and are transmitted via wire 302 toprocessing system 101.

This happens several times a second and allows a picture of the areaimmediately below probe 102 in the patient's body to be built up andshown to the operator on display unit 106.

FIG. 4

Transducer array 303 is shown in more detail in FIG. 4. It comprises aplurality of individual transducers, such as transducers 401, 402, 403,404 and 405. Not all the transducers are shown in this figure for thepurposes of clarity, as indicated by dotted line 406. In this example,probe 102 includes 128 individual transducers.

FIG. 4 also shows that wire 302 is actually made up of 128 separatewires, each controlling one transducer. Thus for example, transducer 401is controlled by wire 411, transducer 402 is controlled by wire 412,transducer 405 is controlled by wire 415, and so on as indicated bydotted line 407. The instructions received from processing system 101are simply instructions on a particular wire for the correspondingtransducer to fire, and thus the probe does not include any processingcapability of its own. The reflection received by each transducer isalso sent along its own wire.

Typically, modern probes use transducer arrays such as that shown inFIG. 4 to produced focussed ultrasound beams. This allows a sweepingalong the array in the lateral direction shown by arrow 408, thusproducing a two-dimensional image by combining the reflections receivedfrom individual transducers. However, other types of probes may usetransducers that are mechanically swept to create a focussed beam, ormay use fewer or more transducers in their arrays, or may usetwo-dimensional arrays. The invention described herein can be used withany type of probe that forms focussed ultrasound beams, as long as thecharacteristics of the probe can be modelled.

FIG. 5

FIG. 5 illustrates how a subset of the transducers in array 303 are usedto create a single ultrasound beam. In this example, a subset 501comprises 32 adjacent individual transducers. These are excited in aparticular order, which will be discussed further with reference to FIG.6, in order to produce a focussed ultrasound beam 502. Typically, theleft-most subset of transducers is fired first, and a step along, onetransducer at a time, is made until the right-most set is fired. Thiscreates a sweep in the direction of arrow 408.

As shown in FIG. 5, beam 502 is focussed at it narrowest point 503. Thisis where the image produced will be sharpest. Prior art methods ofanalysing the information received from probe 102 tend to produceblurring at other depths outside the focal area.

FIG. 6

FIG. 6 illustrates how a subset 501 of transducers produces a singlefocussed beam 502.

The subset comprises, in this example, 32 transducers. Transducer 601 isthe left-most transducer in the set, 602 is the next transducer and soon. Transducer 632 is the right-most transducer in the set withtransducer 631 being adjacent to it. Transducers 616 and 617 are thecentral transducers.

In order to generate the beam 502 a wave front must be generated thatwill converge at the focus 503. This is achieved by first firingtransducers 601, and 632, then transducers 602 and 631, and so on untilfinally transducer 616 is fired with transducer 617. This creates a wavefront having the shape 633 shown in FIG. 6.

In prior art systems, the information returned by transducers 601 to 632would have a delay and sum process applied to it in order to combine asingle time trace for each step along array 303. This results in a smallnumber of time traces that can be easily reproduced as image data andoutput to display 106. However, as previously described, the resultingimage is only sharp at the depth represented by the focus 503, and isblurred everywhere else. The present invention provides a method ofanalysing the data received from all the transducers in array 303 havingundergone a delay and sum process, to provide a less blurred image.

This invention has application for ultrasound probes having differentarrangements of transducers, in addition to that illustrated herein.

FIG. 7

Steps carried out by the ultrasound processing system 101 are describedin FIG. 7. At step 701 the machine is switched on and at step 702instructions are loaded into memory. These could be pre-installedinstructions stored on a hard drive, instructions downloaded from anetworked location, or instructions installed from a local removabledrive, such as a CD-ROM or flash drive.

At step 703 ultrasound images are produced and displayed, and at step704 the machine is switched off.

FIG. 8

Operation of the ultrasound machine shown in FIG. 1 is illustrateddiagrammatically in FIG. 8. The control unit 201 sends pulses 801 to theprobe 102. The probe returns data which is supplied as time traces 802to processing unit 207. Unit 207 supplies image data 803 to display 106and control data 804 to control unit 201, indicating that the cycle maybegin again. It will be appreciated that this is a greatly simplifieddescription of the operation of the machine shown in FIG. 1. Thisdescription will not go into detail regarding the control unit 201, thegeneration of pulses or the receiving of information from the probe. Norwill the beam-forming techniques used to create a beam such as beam 503be further discussed. The invention herein described relates to theprocessing of the time traces 802 and conversion of them into image data803 by processing unit 207. This can be done with any method ofproducing the time traces, by any suitable ultrasound machine and probe.

FIG. 9

FIG. 9 details step 703 at which ultrasound images are produced anddisplayed. At step 901 at least one pulse signal is produced and sent toat least one transducer, and at step 902 answering signals are received.At step 903 these signals are converted to time traces and at step 904 aquestion is asked as to whether the sweep of array 303 is complete. Ifthis question is answered in the negative, control is returned to step901 and more pulses are generated.

However, if the sweep is complete then at step 905 image data isproduced from the time traces by processing unit 207 and at step 906 itis displayed to the operator on display 106.

Control is then returned to step 901. However, at any time an interruptmay be received indicating that the machine is to be switched off,leading to a completion of step 703.

The remainder of this description will concentrate on further detailingstep 905.

FIG. 10

FIG. 10 illustrates the image data being produced by the ultrasoundmachine shown in FIG. 1. It includes an area 1001 corresponding to thearea of the patient being imaged. The top of the picture represents thearea closest to the probe 102.

Area 1002, marked with dotted lines, represents the area in which atraditional ultrasound image would be focussed. The remainder of theimage, constituting most of the image space 1001, would be blurred to agreater or lesser degree. The invention herein described aims to improveupon traditional ultrasound by providing a less blurred imagethroughout.

Axes 1003 illustrate the reference frame of the area being imaged: r₁represents the lateral direction, or the direction along transducerarray 303. Direction r₂ represents the elevation direction and directionr₃ is the depth into the specimen. Any point, such as point 1004, inspace 1001 can therefore be represented as a three-dimensional point(r₁, r₂, r₃). When an ultrasound beam from a transducer is reflectedfrom it, this generates a time trace, which is some function of r₁, r₂and t, where t is the time taken to receive an echo at the transducer.

FIG. 11

FIG. 11 is a representation of how the reflection of a point in space isconverted to a time trace by a transducer. The point 1101 reflects anultrasound pulse 1102 generated by transducer 1103. This results in atime trace 1104 which has co-ordinates r₁ and r₂ in the lateral andelevational directions respectively, and a first time measurement t₁,representing the amount of time it took for the pulse to be reflected. Apoint 1105 which is on the same line as point 1101 but at a smallerdepth, would have the same lateral and elevational entries in its timetrace 1106, but a smaller time measurement t₂. In a perfect world, itwould therefore be a simple matter to equate time with depth and displaythe time traces to an operator accordingly. However, a pulse is not aneat line, but a wave having a focal point. In reality, each scattereris reflected back to more than one transducer.

FIG. 12

FIG. 12 is a more real-world illustration of what happens when a pulseis reflected back from a scatterer. Transducer 1103 produces a pulse1201 that has a waveform having a focus at point 1202. When a pluralityof transducers are placed side-by-side, each produces a pulse thatinterferes with its neighbours' waveforms.

The pulses of transducers 1203 and 1204, which are adjacent totransducer 1103, are also shown in dotted lines. The delay and sumtechnique commonly used by ultrasound probes takes advantage of theHuygens principle, whereby each transducer's pulses are considered as anindividual source that contribute to an overall wave front. Asillustrated in the Figure, the beams produced by transducers 1203, 1103and 1104 therefore contribute to an overall wave front having a focalpoint defined where the waves sum coherently. Thus, waves reflected by ascatterer located at point 1101 will also be coherent, and when eachtransducer's time trace is displayed, the image of the scatterer willappear sharp.

However, point 1105 is closer to the transducers and so is outside thefocal area. The time trace produced by each transducer is thereforeslightly different in both the lateral and the time directions, asreflections from a scatterer at point 1105 are not coherent. Thus, asimple reproduction of the time traces as image data will lead to ablurring at point 1102, i.e. it will look larger than it should be.

When this effect is taken into account for the many scatterers that arepresent in a specimen being imaged, then it will be understood that ifmany scatterers are enlarged through blurring then they will start tomerge together, and a loss of sharpness outside the focus depth is theresult.

FIG. 13

The described invention therefore proposes a different way of dealingwith the time traces, rather than traditional delay and sumpost-processing techniques which fail to take account of the blurring.

For any given probe and for a type of specimen being imaged, it ispossible to construct a blurring matrix which encapsulates how ascatterer at each point in the space being imaged will be reflected backto the transducers in the array. This is a function of a number ofparameters, such as the way in which the transducers are fired, thespeed of sound through the specimen being imaged, and so on. Thus, foreach probe and type of specimen, a blurring matrix can be calculated orretrieved from memory as necessary. Thus the ultrasound process can berepresented as shown in FIG. 13.

The field of scatterers can be modelled as a vector 1301, x. This ispassed through a blurring matrix 1302, H, and is also subject to noise1303, n, to produce time traces 1304, y. This leads to equation 1310,that y is equal to the sum of the product of H and x, and n:

y=Hx+n

Given that y is known, H is known and n can be modelled, we can arriveat an estimate for x. This estimate is what is displayed on display 106as image data, and is as close as possible to the original scattererfield.

However, estimating x is not simple because H is very large and is notinvertible. Typically H, a non-sparse matrix, might have 1,000,000 rowsand columns. Thus, obtaining x using a simple multiplication by theinverse of H is not possible. Further, any calculation using H will beslow and may not fit in the memory of the processing unit.

Thus the problem is to estimate x given y and H. In order to do this, xis modelled as a product of its echogenicity 1306, S, and a randomnoise-like component 1307, w:

x=Sw

S is a diagonal matrix of real valued, non-negative echogenicity valuesand w is a random vector of complex valued zero mean Gaussian values. Sis the same size as H and w is the same size as x. However, because S isdiagonal it is easier to use than H.

The noise 1303 is assumed to have a circularly symmetric complexGaussian distribution with zero mean and co-variance Σ_(n), as shown inequation 1308:

n˜CN(0,Σ_(n))

Given this, the estimation of x can be configured as inferring thedistribution of x conditioned on the data y, but with consideration tothe parameters S and Σ_(n). This can be done by iterating betweenfinding the distribution of x and improving point estimates of S andΣ_(n). The best estimate 1309 of x can be written as follows:

{circumflex over (x)}=(H ^(H)Σ_(n) ⁻¹ H+S ⁻²)⁻¹ H ^(H)Σ_(n) ⁻¹ y

This is the maximum a posteriori estimate for x, which is equivalent tothe mean of the conditional distribution on x.

FIG. 14

By making certain assumptions about the underlying scatterer field,equation 1309 can be simplified and made more calculable.

First, the nature of the noise 1303 is considered. Noise may come from avariety of sources. For example, white noise is introduced at theamplifier stage, cross channel coupling between transducer wires canoccur, and so on. Because so much of the noise is electrical, prior artsystems generally assume that the noise is invariant across thescatterer field, i.e. the time trace from each transducer includes thesame amount of noise.

However, noise is also due to other factors. Non-linear effects can becaused by pulses being reflected off more than one scatterer, or probeinterface reverberation. Noise may also be caused by reflection to ascatter in the side lobes on an ultrasound beam or simply an incorrectcalculation of the blurring matrix 1302. Consideration of this type ofnoise gives rise to the understanding that the noise field varies acrossthe scatterer field, and may in fact be scaled across the imagedependent upon the scatterer field.

Thus an assumption can be made that the co-variance Σ_(n) of the noiseis in proportion to the co-variance S² of the non-conditionaldistribution of x. This can be written as equation 1401, that theproduct of Σ_(n) and the inverse of S² is some constant q multiplied bythe identity matrix:

Σ_(n) S ⁻² ηI

A value of 0.3 for η has been found to be effective.

While the assumption is made that the noise variance varies with thesignal variance, another assumption can be made that the noise variancevaries slowly, such that it can reasonably be assumed to be constantacross the support of the point spread function. This leads to theapproximation that the product of the Hermitian transpose of H and theinverse of Σ_(n) is approximately equal to the product of the inverse ofΣ_(n) and the Hermitian transpose of H, as shown at 1402:

H ^(H)Σ_(n) ⁻¹≈Σ_(n) ⁻¹ H ^(H)

As described with respect to FIG. 13, a technique for finding anestimate of the image data x is to iterate to a solution of the mean ofthe conditional distribution on x; the kth iteration of this is writtenas m_(k). The kth estimate m_(k), following equation 1309, is thereforethe product of the following: the inverse of the sum of the product ofthe Hermitian transpose of H, the inverse of the kth iteration of Σ_(n)and H, and the kth iteration of the inverse of S²; the Hermitiantranspose of H; the kth iteration of the inverse of Σ_(n); and y, asshown at 1403:

m _(k)=(H ^(H)Σ_(n,k) ⁻ H+S _(k) ⁻²)⁻¹ H ^(H)Σ_(n,k) ⁻¹ y

Substituting in approximation 1402 gives the simplified estimate m_(k)to be equal to the product of the following: the inverse of the productof the Hermitian transpose of H and H, added to the product of the kthiteration of Σ_(n) and the kth iteration of the inverse of S²; theHermitian transpose of H; and y, as shown at 1404:

m _(k)=(H ^(H) H+Σ _(n,k) S _(k) ⁻²)⁻¹ H ^(H) y

Using assumption 1401, an initial estimate for m can therefore bewritten as the product of the following: the inverse of the product ofthe Hermitian transpose of H and H, added to the product of η and theidentity matrix; the Hermitian transpose of H; and y, as shown at 1405:

m ₀=(H ^(H) H+ηI)⁻¹ H ^(H) y

Thus an initial approximation m₀ can be made that is dependent only uponthe blurring matrix 1302, H, and the time traces 1304, y. This is stilla very large calculation, but, as will be shown in FIG. 15, if it isconsidered in the appropriate domain it becomes tractable.

FIG. 15

As previously discussed, both Σ_(n) and S are diagonal matrices, andthus, although large, are sparse and can fit in memory. However, theblurring matrix H is very large and full, and calculations with itrequire a large amount of memory and time.

So far only the lateral time domain 1501 has been considered, which hasaxis in the lateral r₁ direction, the elevation r₂ direction and time. AFourier transform 1502 in the r₁ and r₂ directions transforms domain1501 into the lateral K-time domain 1503, where the lateral direction istransformed into the K-lateral direction and the elevation directioninto the K-elevation direction respectively. Because the blurring matrix1302, H, is invariant across the lateral direction in domain 1501,applying Fourier transform 1502 to it gives a block diagonal matrix1504. A similar effect occurs with the elevational direction. Matrix1504 is made up of a plurality of sub-matrices, such as matrices 1508,1509, and so on. Each of these sub-matrices is small and can thereforefit into the memory of a typical processing system.

Thus, in order to calculate, for example, matrix 1505, which is thefirst part of equation 1405, a plurality of sub-matrices 1506 can becalculated and combined into a single matrix, the inverse Fouriertransform of which is the required matrix. In other words, in order tocalculate the inverse of the sum of the product of the Hermitiantranspose of H and H and the identity matrix scaled by q:

(H ^(H) H+ηI)⁻¹

the Fourier transform 1502 is applied to the matrix H to produce aplurality of sub-matrices H_(sub,i), where i indicates the ithsub-matrix. Because ηI is a diagonal constant matrix, it passesunchanged through the Fourier transform 1502. Thus, for each value of i,the inverse of the sum of the product of the Hermitian transpose ofH_(sub,i) and H_(sub,i), and the appropriate sub-matrix of the identitymatrix scaled by η, can be calculated:

(H _(sub,i) ^(H) H _(sub,i) +ηI _(sub,i)) ⁻¹

The resultant sub-matrices are combined into a single matrix andtransformed back to the lateral time domain 1501 using inverse Fouriertransform 1507 to give a result for matrix 1505.

Other calculations involving the blurring matrix H can be similarlycarried out in domain 1503. However, because Σ_(n) and S are diagonal inthe lateral time domain, calculations involving these matrices arecarried out in domain 1501.

FIG. 16

The procedure for estimating the scatterer field x is showndiagrammatically in FIG. 16.

A first estimate 1601 is made of the distribution of the image data.This estimate is used to make a first estimate 1602 of the echogenicity,S₀, and a first estimate 1603 of the co-variance of the noise, Σ_(n,0).The method of obtaining these last two estimates will be describedfurther with reference to FIGS. 20 to 22.

The first estimates of the echogenicity and of the variance of the noiseare then used to update the estimate of the mean of the distribution ofthe image data at 1604. This estimate is then used to update theestimates for the echogenicity and the co-variance of the noise at 1605.The process then iterates between step 1604 and step 1605 until asatisfactory convergence is reached at 1606. The threshold forconvergence can be defined as necessary. However, around six iterationsof steps 1604 and 1605 are generally required to reach a satisfactoryconvergence. This algorithm is a form of the Expectation Maximisationalgorithm.

The first estimate of the mean of the distribution of the image data m₀found at 1601, is calculated using equation 1405. This involvescalculating matrix 1505 as shown in FIG. 15. This matrix is stored forlater use as it is used in the step of updating the estimates m_(k),using equation 1403, at step 1604.

In an alternative embodiment, first estimates of the echogenicity andthe co-variance of the noise can be obtained first and used in step 1604to obtain a first estimate of m_(k).

In this description, the input data into the above procedure is the timetraces, y, and the ultimate result will be dependent upon this inputdata. However, it is possible that the procedure could be used tode-blur other, similar, data, and thus the input data is not limited totime traces.

FIG. 17

FIG. 17 details step 905 carried out by processing system 101 to produceimage data from the input time traces. This uses the process outlined inFIG. 16.

First at step 1701 an appropriate blurring matrix H and constant q areretrieved from memory. The variable k is set to zero. At step 1702 thematrix 1505 is calculated, and at step 1703 this matrix, H and the inputtime traces are used in equation 1405 to generate the first estimate m₀of the image data.

At step 1704 the estimate m₀ is used to calculate the first estimate S₀of the echogenicity, using any appropriate method. At step 1705 it isused to calculate the first estimate Σ_(n,0) of the co-variance of thenoise, as will be further described with reference to FIG. 20.

At step 1706 k is incremented by one and, at step 1707, what are nowestimates S_(k−1) and Σ_(n,k−1) are used to update the estimate m_(k) ofthe image data using equation 1404, At step 1708 this new estimate m_(k)is used to update the estimate S_(k) of the echogenicity and Σ_(n,k) ofthe covariance of the noise. At 1709 a question is asked as to whetherthe convergence of the process is satisfactory. If this question isanswered in the negative then control is returned to step 1706, k isincremented by one and steps 1707 and 1708 are carried out again. If,however, the convergence is satisfactory then step 905 is completed. Thelast estimate m_(k) of the mean of the distribution of the scattererfield, is used as the best estimate of the scatterer field, and istherefore output as image data to be displayed on display unit 106. Themethods of performing the calculations of m_(k), S_(k) and Σ_(n,k) atsteps 1707 and 1708 will now be detailed with respect to FIGS. 17 to 22.

FIG. 18

As described with reference to FIG. 15, the initial estimate of theimage data m₀ can be calculated in the lateral K-space time domain 1503,because it includes only the blurring matrix H, which is diagonalised bythe Fourier transform, and a constant. Matrix 1403, which includes theterms Σ_(n) and S, cannot be calculated in this way, because Σ_(n) andS, while diagonal in the lateral time domain 1501, become non-diagonaland full in the lateral K-space time domain 1503. Thus, a differentmethod must be used to calculate m_(k).

The Conjugate Gradients algorithm is used at each iteration of step1707. Equation 1403 can be rewritten as equation 1801, as follows:

${\overset{\overset{A}{}}{\left( {{H^{H}{\sum\limits_{n,k}^{- 1}H}} + S_{k}^{- 2}} \right)}\overset{\overset{c}{}}{m_{k}}} = \overset{\overset{b}{}}{H^{H}{\sum\limits_{n,k}^{- 1}y}}$

which is of the form Ac=b, as shown at 1802. The Conjugate Gradientsalgorithm solves this type of equation by finding the value of c thatminimises the least squares difference. It is guaranteed to converge inthe same number of iterations as the dimension of the matrix A. However,in our case the matrix A may have dimensions of order 1,000,000 by1,000,000 entries, which makes for very slow convergence.

The algorithm can be speeded up by use of a preconditioner P. Theequation to be solved is rewritten as AP−1 Pc=b, as shown at 1803. IfP−1 is cheap to compute, the Conjugate Gradients algorithm can then beused to solve for Pc, from which c can be found. The idea ofpreconditioning is to produce the matrix AP−1 that is more amenable tofast inversion than A on its own, meaning that the algorithm convergesin fewer iterations. A good preconditioner is an approximation to A−1that is cheaply invertible.

In the case under consideration, the matrix A−1 is the first part ofequation 1403, which as discussed with reference to FIGS. 14 and 15 canbe approximated by matrix 1505. At this stage of the process, thismatrix has already been calculated in order to provide the firstestimate of the image data. Not only is it a close approximation of therequired matrix, but it is relatively easy to invert. Thus it can beused as the preconditioner in the Conjugate Gradients algorithm, asshown at 1804.

FIG. 19

FIG. 19 details step 1707, at which the next estimate m_(k) of the imagedata is obtained using the conjugate gradients algorithm described withrespect to FIG. 18.

At step 1901 the matrices P⁻¹ (matrix 1505), c₀ (the current estimatem_(k) of the image data) and A (matrix 1406) are retrieved from memoryand at step 1902 a number of variables are initialised as follows: r₀ isinitialised to be the product of A and c₀ subtracted from b; z₀ isinitialised to be the product of P⁻¹ and r₀, p₀ is initialised to be z₀,and a counter n is set to zero:

r ₀ =b−Ac ₀

z ₀ =P ⁻¹ r ₀

p ₀ =z ₀

n=o

At step 1903 the variables are updated as follows: a variable α is setto be the product of z_(n) and r_(n) divided by the product of theHermitian transpose of p_(n), the Hermitian transpose of A and p_(n);c_(n+1) is updated to the sum of c, and the product of α, A and p_(n);r_(n+1) is updated to the product of α, A and p_(n) subtracted fromr_(n); a variable β is set to be the product of the Hermitian transposeof z_(n+1) and r_(n+1) divided by the product of the Hermitian transposeof z_(n) and r_(n); and p_(n+1), is updated to be sum of z_(n+1) and theproduct of β and p_(n):

$\alpha = \frac{z_{n}r_{n}}{p_{n}^{H}A^{H}p_{n}}$c_(n + 1) = c_(n) + α p_(n) r_(n + 1) = r_(n) − α Ap_(n)z_(n + 1) = P⁻¹r_(n + 1)$\beta = \frac{z_{n + 1}^{H}r_{n + 1}}{z_{n}^{H}r_{n}}$p_(n + 1) = z_(n + 1) + β p_(n)

At step 1904 n is incremented by one, and at step 1905 a question isasked as to whether satisfactory convergence has been reached, indicatedby the closeness of the calculated r_(n+1) to zero. If it is not closeenough, another iteration of step 1903 is carried out. Alternatively,convergence is satisfactory and step 1707 is complete. The last value ofcn+1 is used as the result and becomes the next estimate of mk.

Other methods of implementing the Conjugate Gradients algorithm areavailable, but in all cases the pre-conditioner 1804 may be used tospeed up the convergence.

FIG. 20

The description now turns to the way in which the estimate of thecovariance of the noise Σn is initialised and updated during steps 1705and 1708 respectively.

As described with reference to FIG. 14, prior art methods of analysingultrasound data assume that the noise 1303, n, is invariant across thescatterer field, but better results can be obtained by assuming thatthis is not the case. However, around a scatterer it is assumed thatthere exists some region where the noise is effectively constant. Thus,for example, in our image space 1001 there is a region 2001 around ascatterer 2002 in which the noise can be assumed to be substantiallyconstant. Starting from equation 1310, that the time trace is equal tothe product of the blurring matrix and the image data minus the noise,the noise can be written as the product of the blurring matrix and theimage data subtracted from the time trace, as shown at 2003.

Thus the variance of the noise within region 2001 can be estimated usingthe sum of the squares of the residuals within that region, using thelatest estimate m_(k) of the image data. This can be done for aplurality of regions to obtain a result for the estimate of thecovariance matrix Σ_(n) at the kth iteration.

The procedure is simplified by assuming that Σ_(n,k) is diagonal,meaning that only the diagonal elements need to be calculated. Theelement in the ith row and ith column is estimated as follows. A set ofentries in the current estimate of the image data, m_(k), is denotedS_(sub,i). This corresponds to a region such as region 2001 of the imagedata. For each point in this region, the magnitude of the residual ofthe blurring matrix H multiplied by the current estimate m_(k) of theimage data, subtracted from the time trace y is squared. The result isdivided by two and summed for all j in the set S_(sub,i), and this sumis divided by the size of the set S_(sub,i). This is given by equation2004:

$\left( {\hat{\sum}}_{n,k} \right)_{i,i} = {\frac{1}{N_{S_{{sub},i}}}{\sum\limits_{j \in S_{{sub},i}}{\frac{1}{2}{\left( {y - {Hm}_{k}} \right)_{j}}^{2}}}}$

Each entry in the diagonal matrix Σ_(n,k) is calculated like this,giving the full estimate Σ_(n,k) at step 1708. This estimate of thevariance of the noise is made possible by assuming that the noise variesacross the image data, but is locally invariant in any small enoughregion.

This equation is used to initialise the estimate of the variance of thenoise using m0, as well as to update it at each iteration.

This method of estimating the variance of the noise could be used withother versions of the Expectation Maximisation algorithm. In particular,it could be used with other methods for initialising and updating boththe estimate of the mean of the distribution of the image data and theestimate of the echogenicity.

The echogenicity S_(k) is initialised and updated at each iteration inany suitable way.

FIG. 21

The algorithm described herein with reference to FIGS. 13 to 20estimates image data to be displayed on a display, given input data as aplurality of time traces and a blurring matrix. It can be used toprovide less blurred images than those currently available usingtraditional ultrasound machines. However, it is also useful in that itcan provide good quality ultrasound images using a relatively smallamount of processing power memory.

A portable ultrasound machine is illustrated in FIG. 21. Portablemachine 2101 resembles a laptop computer, and has similar limitations interms of the available processing power and memory. Typically, a probe2102 plugs into a port 2103, such as a USB port, and the time traces itproduces are displayed on display 2104. A keyboard interface 2105 isprovided to allow control of the system by an operator.

In an embodiment of the invention, a device 2106 is placed between probe2102 and port 2103. The USB device configures system 2101 to carry outthe method described herein, rather than its normal method. Thus, astandard portable machine can be easily modified to produce improvedimage data. Other methods of introducing executable instructions intothe processor of an ultrasound machine are also contemplated.Alternatively, newer ultrasound machines may as standard be built withthese instructions embedded.

What is claimed is:
 1. Apparatus for producing image data representing aspecimen being imaged using ultrasound, comprising a processor andmemory, wherein said processor is configured to: receive input datarepresenting an output from an ultrasound probe; process said input datato produce output image data; and output said output image data; saidprocessor being configured to carry out said step of processing saidinput data by: obtaining a first matrix that is dependent uponattributes of said probe and on the way in which said input data wasproduced; initialising estimates for the image data, the echogenicity ofsaid image data, and the variance of the noise of said image data;iteratively performing the following steps until a predetermined levelof convergence is reached: (a) update said estimate of said image datausing: said first matrix, a respective said estimate of saidechogenicity, a respective said estimate of said variance of said noise,and said input data, (b) update said estimate of the echogenicity usingsaid updated estimate of said image data, and (c) update said estimateof the variance of the noise using said updated estimate of said imagedata; wherein step (c) is performed using an update rule that assumesthe noise to be varying across said image data but locally invariant forany small enough region of said image data; and using a convergedestimate of said image data to produce said output image data. 2.Apparatus according to claim 1, wherein said update rule for step (c)is:$\left( {\hat{\sum}}_{n,k} \right)_{i,i} = {\frac{1}{N_{S_{{sub},i}}}{\sum\limits_{j \in S_{{sub},i}}{\frac{1}{2}{\left( {y - {Hm}_{k}} \right)_{j}}^{2}}}}$where: ({circumflex over (Σ)}_(n,k))_(i,i) is an entry in the ith rowand ith column of a matrix representing the estimate of the variance ofthe noise at a kth iteration; H is said first matrix; m_(k) is theestimate of the image data at the kth iteration; y is said input data;S_(sub,i) is a set of entries in m_(k) corresponding to a region of saidimage data across which the noise is assumed to be invariant; and N_(S)_(sub,i) is the size of S_(sub,i).
 3. Apparatus according to claim 1,wherein said processor is configured to initialise said estimate of theimage data using said first matrix and said image data.
 4. Apparatusaccording to claim 3, wherein said processor is further configured toobtain a constant and evaluate a second matrix using said first matrixand said constant, and is configured to initialise said estimate of theimage data additionally using said second matrix.
 5. Apparatus accordingto claim 4, wherein said processor is configured to initialise saidestimate of said image data by calculating the product of: (a) saidsecond matrix, (b) a transpose of said first matrix, and (c) said imagedata.
 6. Apparatus according to claim 4, wherein said second matrix isthe inverse of the sum of: (a) the product of: a transpose of the firstmatrix, and the first matrix, and (b) the product of: said constant, andan identity matrix having the same dimensions as said first matrix. 7.Apparatus according to claim 4, wherein said constant representsnoise-to-signal ratio of the image data.
 8. Apparatus according to claim4, wherein an algorithm used in step (a) is the Conjugate Gradientsalgorithm.
 9. Apparatus according to claim 8, wherein each iteration ofstep (a) uses said second matrix as a preconditioner in the algorithm.10. Apparatus according to claim 1, wherein said processor is configuredto update the estimate of the echogenicity additionally using said inputdata.
 11. A method of producing image data representing a specimen beingimaged using ultrasound, comprising the steps of: receiving input datarepresenting an output from an ultrasound probe; processing said inputdata to produce output image data; and outputting said output image datato a display; wherein said step of processing said input data comprises:obtaining a first matrix that is dependent upon attributes of said probeand on the way in which said input data was produced; initialisingestimates for the image data, echogenicity of said image data, and thevariance of noise of said image data; iteratively performing thefollowing steps until a predetermined level of convergence is reached:(a) update said estimate of said image data using: said first matrix, arespective said estimate of said echogenicity, a respective saidestimate of said variance of said noise, and said input data, (b) updatesaid estimate of the echogenicity using said estimate of said imagedata, and (c) update said estimate of the variance of the noise usingsaid estimate of said image data; wherein step (c) is performed using anupdate rule that assumes the noise to be varying across said image databut locally invariant for any small enough region of said image data;and using a converged estimate of said image data to produce said outputimage data.
 12. A method according to claim 11, wherein said update rulefor step (c) is:$\left( {\hat{\sum}}_{n,k} \right)_{i,i} = {\frac{1}{N_{S_{{sub},i}}}{\sum\limits_{j \in S_{{sub},i}}{\frac{1}{2}{\left( {y - {Hm}_{k}} \right)_{j}}^{2}}}}$where: ({circumflex over (Σ)}_(n,k))_(i,i) is an entry in the ith rowand ith column of matrix representing the estimate of the variance ofthe noise at a kth iteration; H is said first matrix; m_(k) is theestimate of the image data at the kth iteration; y is said input data;S_(sub,i) is a set of entries in m_(k) corresponding to a region of saidimage data across which the noise is assumed to be invariant; and N_(S)_(sub,i) is the size of S_(sub,i).
 13. A method according to claim 11,wherein said step of initialising said estimate of the image data iscarried out using said first matrix and said image data.
 14. A methodaccording to claim 13, further comprising the steps of obtaining aconstant and evaluating a second matrix using said first matrix and saidconstant, wherein said step of initialising said estimate of the imagedata is carried out additionally using said second matrix.
 15. A methodaccording to claim 14, wherein said step of initialising said estimateof said image data is carried out by calculating the product of: (a)said second matrix, (b) a transpose of said first matrix, and (c) saidimage data.
 16. A method according to claim 14, wherein said secondmatrix is the inverse of the sum of: (a) the product of a transpose ofthe first matrix, and the first matrix, and (b) the product of saidconstant and an identity matrix having the same dimensions as said firstmatrix.
 17. A method according to claim 14, wherein said constantrepresents noise-to-signal ratio of the image data.
 18. A methodaccording to claim 14, wherein an algorithm used in step (a) is theConjugate Gradients algorithm.
 19. A method according to claim 18,wherein each iteration of step (a) uses said second matrix as apreconditioner in the algorithm.
 20. A non-transitory computer-readablemedium encoded with program instructions executable by a computer that,when executed by a computer, cause a computer to perform steps of:receiving input data representing an output from an ultrasound probe;processing said input data to produce output image data; and outputtingsaid output image data to a display; wherein said step of processingsaid input data comprises: obtaining a first matrix that is dependentupon attributes of said probe and on the way in which said input datawas produced; initialising estimates for the image data, theechogenicity of said image data, and the variance of the noise of saidimage data; iteratively performing the following steps until apredetermined level of convergence is reached: (a) update said estimateof said image data using: said first matrix, a respective said estimateof said echogenicity, a respective said estimate of said variance ofsaid noise, and said input data, (b) update said estimate of theechogenicity using said estimate of said image data, and (c) update saidestimate of the variance of the noise using said estimate of said imagedata; wherein step (c) is performed using an update rule that assumesthe noise to be varying across said image data but locally invariant forany small enough region of said image data; and using said convergedestimate of said image data to produce said output image data.